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A practical method for computing the effectiveness of a space nuclear shield perforated by small 

tubing and cavities has been developed and tested. 

Analysis of ducts and voids in a shield cannot be performed using the two-dimensional transport 
methods that are presently the most efficient means of analyzing reactor and shield proolems. A three- 
dimensional Monte Carlo analysis is the only rigorous approach; however, it is too costly to consider in 
tradeoff studies. Methods are needed that will provide reasonable accuracy at a moderate cost. 

Our approach to the problem is to use solutions of a two-dimensional transport code and eval- 
uate perturbations of that solution using last flight estimates and other kernel integration techniques. In 
general, perturbations are viewed as a change in the source strength of scattered radiation and a change 
in attenuation properties of the region, and the effect of these changes is evaluated by performing last 
flight estimates and point kernel integrations. This method has been implemented in a new computer 
program called DUCT. 

Calculations performed with the DUCT code have been compared to experiments conducted at the 
Gulf Radiation Technology Linac facility. The shield studies involved measurements on a lithium hy- 
dride assembly with a sodium filled duct and a NERVA-BATH shield mockup. Comparison has also 
been made with the MORSE Monte Carlo code. Overall agreement between measurements, DUCT, and 
MORSE has been better than 10 percent. 


1. INTRODUCTION 

A variety of techniques are available for 
estimating radiation streaming down ducts. Some 
of these techniques are summarized in ref. 1. 

The majority of. the duct streaming work, both : 
experimental and analytical, has been concentrat- 
ed on void ducts whose dimensions a re com- 

pared to the radiation mean free path in the mat- 
erial. For these cases, the radiation incident on 
the mouth of the duct is the predominant source 
of radiation reaching the duct exit either by di- 
rect streaming or by albedo from the duct wall. 
Relatively simplified analyses of these two com- 
ponents provide answers that are in reasonable 
agreement with experiment for cases such as 
shelter entry ways or beam ports through reactor 
shields. 

Ducts through shields in space reactors 
are, however, more likely to be filled either by 
metal or by some fluid so that radiative exchange 
between the duct material and the shield material 
becomes an important contribution to the radia- 
tion exiting from the duct. Evaluating the effect 
of a filled duct is therefore a more complicated 
process and is not amenable to an unsophisticated 
approach. An exact solution requires the use of 
a radiation transport code with a three-dimension- 
al generalized geometry capability (a Monte Carlo 
code). The motive behind the work described 
here is to arrive at an approximate method that 
can be applied to iterative design studies without 


the large expenditure of computer time required 
for a series of Monte Carlo calculations. 

This research was divided into four parts; 
namely, 

1. Formulation and programming of a code 
for calculating the effects of shield irre- 
gularities on the space craft environ- 
ment. This method consists of using the 
solution of a two-dimensional transport 
code as a starting point and making per- 
turbations to that solution to account for 
shield irregularities such as coolant ducts 
and other shield perturbations. 

2. Performance of analytical tests to deter- 
mine the range of validity of the code 
developed in Part 1. 

3. Performance of an experimental program 
using a geometry and material similar to 
real space-power shielding problems but 
'’clean" enough to provide a good check on 
the analysis. 

4. Extension of the experimental-analytical 
comparison to include neutron and second- 
ary gamma production and transport in a 
NERVA-BATH shield assembly. 

2. DEVELOPMENT OF THE DUCT CODE 

Except for isolated irregularities, we 
assume that the shield can be approximated by 
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use of a two-dimensional model. Neglecting 
these irregularities, a two-dimensional calcula- 
tion can be performed using a transport code such 
as TWOTRAN (ref. 2) or DOT (ref. 3). These 
codes can be made to generate a M flux tape" 
which is a complete listing of fluxes within a spa- 
tial energy and angle mesh specified in the calcu- 
lation. The DUCT code that we have written takes 
these fluxes as input and solves for perturbations 
in the flux at specified detector points caused by 
the presence of irregularities in the shield. In the 
present form, the code only considers cylindrical 
ducts. In some parts of the code the 2D angular 
flux mesh is used in a discrete manner, in other 
parts, interpolations in angle are made. Spatial 
interpolations are made in all cases and the fixed 
energy grouping is used throughout. 

The primary assumption attendant to the 
derivation of the method is that the inward- di- 
rected current passing through the surface of a 
duct is obtained by integrating over the angular 
flux predicted from the unducted 2D solution. In 
other words, it is assumed that changing the mat- 
erial in the duct region does not change the level 
of radiation passing into the region. Auxilliary 
one-dimensional transport calculations in cylin- 
drical geometry show this to be a good approxi- 
mation. 

2.1 ESTIMATOR FOR THE SCATTERED 
RADIATION 


The last collision kernel volume integral 
is given as 


9<E) S = 
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where 

cp(Q' } r, 2, E 1 ) are flux components in phase space 
coordinates in cylindrical geometry cen- 
tered in dv. 

O’ is direction of incoming flux, and E' is en- 
ergy of incoming flux. 

R is the distance from dv to the detector. 

0 is the direction from dv to the detector. 

E is the energy after collision. 

P(E'0 Eft) is the scattering cross section per 
steradian per cc for transitions from O' 
to and E 1 to E. These two transitions are 
not independent since for a given scatter- 
ing species the energy lost in a collision 
depends on the scattering angle. Where dv 
is filled with a variety of scattering species 
more than one E can result from a given 
incident E' and Q 1 flux component. 

b(E’) is E t (E , )R. 

To simplify the evaluation of the last 
collision integral, group-averaged cross sections 
and group fluxes are used for both the incident 
and scattered radiation. This allows us to use the 
same cross-section set in the perturbation calcu- 
lation as that used in the 2D transport analysis of 
the idealized shield. 

Another simplification is to restrict the 


the duct shape to cylinders and divide the duct 
volume by planes perpendicular to the duct axis. 
Each of the volume elements can then be treated 
as an equivalent disk source. The uncollided flux 
from an isotropic disk source of strength per 
unit area at a detector on the axis is given by 

S A r 1 

<P=— [^(^t) - E^/itsec 0)J , 

where fut is the integral of total cross sections 
along the axis between the two points and 0 is the 
angle between a ray on the axis and a ray from the 
detector to the outer rim of the disk and the func- 
tion Ej^(x) is given by 
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Using the group averaging approach the flux at a 
detector on the duct centerline at the exit plane is 
then given by 


cp S (g ' ) dCZX} P s (g-g ' >n:,) Vm (a,g) 

g m n 
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where Ps(g-* g'. rn) is the unit volume probability of 
scattering from group go to group g' from the mth 
angle to the direction of the detector. P s (g-*g'm) 
is given by 

n 

p s (g-g ; e g) = ’ 


where 

Sgtg^g 1 ) are given by the group average cross 
sections generated by GGC4, 
is the quadrature weight of the mth angle 
increment = solid angle/4n. 

tp m (n, g) is the angular flux at the nth duct point 
in the mth angle in group g (4ttx neutrons 
cm' 2 STER" 1 ). 

/i(g)R(n) is the total number of mean free paths 
from the nth duct point to the detector for 
group g. 


2. 2 ESTIMATOR FOR THE UNCOLLIDED 
FLUX TO ON- AXIS DETECTOR 

The uncollided flux from radiation passing 
through the duct wall is given by 


9(E) d = 


/*<0, r, Z, E)Q • n dS exp(-E^(E)R) 


where 

cp(Q, r, Z, E) are the flux components. 

0 is the direction from the surface element dS 
to the detector. 

n is a unit vector normal to the surface incre- 
ment. 
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r, Z are the space coordinates centered in dS. 

R is the distance from r, Z to the detector. 

Using the same numerical approach to 
evaluate this integral as was used in the scatter- 
ing calculations detailed above, 


cd (n, g)AS cos Q exp(-/i( g)R(n)) 


where 

AS is the duct surface area associated with the 
particular duct segment. 
cp(n, g) / 4tt is the flux per steradian in the detec- 
tor direction averaged over AS. 

0 is the angle between the surface normal and 
a ray to the detector. 

The perturbation caused by the material 
filling the duct is assumed to express itself as a 
change in the scattering intensity within that geo- 
metric region and as a change in the component 
which streams from the duct wall to the detector 
without suffering a collision. The scattered flux 
is thus obtained by integrating a last flight esti- 
mator over the volume based on a collision den- 
sity given by the unperturbed flux times the dif- 
ferential angular scattering probability. The 
direct flux is obtained by integrating an uncollided 
kernel over the surface of the duct where the local 
source is given by interpolating the ZD fluxes. 

The perturbation is then obtained by comparing 
the results of these calculations with calculations 
where the volume of concern is filled with shield 
material. 

2. 3 ESTIMATOR FOR THE UNCOLLIDED 
COMPONENT FOR OFF-AXIS 
DETECTORS 

Formulations have been made for calculat- 
ing the uncollided component to detectors off the 
duct axis. It will be recalled that the basic me- 
thod is capable of calculating both a collided and 
uncollided component on the duct axis, but for 
off-axis detectors it can only calculate the collid- 
ed contribution. The formulation presented here 
permits an estimate to be made of the enhanced 
penetration of flux vectors whose paths lie parti- 
ally within the duct and partially within the sur- 
rounding shield. 

The complicated geometry is one of the 
major obstacles in formulating the estimator. The 
effective area on the duct wall that can radiate 
toward the detector is given by integrating over 
all the surface areas associated with the duct 
increment whose normal has a component in the 
detector direction 

/ r * n dA , 

AA 

where r is a unit vector in the direction of the 
detector and n is a unit vector normal to the 
tangent plane to the surface. It is assumed that 


the direction 7 is essentially the same for all 
points on the duct increment or that the distance to 
the detector is large when compared to the radius 
of the duct. 

It is possible to evaluate the integral over 
the appropriate area and multiply by the flux vec- 
tor to get a total "source" in the detector direction 
but this is inappropriate because it does not account 
for the varying path lengths through the duct. To 
account properly for the varying path lengths the 
surface must be split into azimuthal increments. 
Some trial and error cases showed that splitting 
the 180° band into six equal 30° increments pro- 
vides a reasonable approximation of the projected 
area. The actual area of each increment is At 
and the flux contribution to the detector is 

-b /R 2 

— — — TTr 1 

cp(E) Q = co(r, E)r • n — e 

where At is the length increment, R is thedetector- 
to-area increment separation distance, and b^ is 
the integral of total cross section along R. 

3. CALCULATIONS OF THE LITHIUM 

HYDRIDE SHIELD EXPERIMENT 

The first test for the DUCT code was an 
attempt to simulate a neutron shield applicable to 
to the SNAP type shield. This shield consists of a 
lithium hydride shield with sodium coolant holes. 

The shield was represented in the DOT cal- 
culation as a cylinder with radius 25. 4 cm in rad- 
ius, 25.4 cm high, having 22 radial intervals and 
22 axial intervals. The top two and bottom two 
axial intervals and the outermost two radial inter- 
vals used steel cross sections; the others were 
LiH. An asymmetric 100-angle quadrature set 
with 76 downward and 24 upward angles was se- 
lected in an attempt to prevent ray effects. The 
spacial mesh was nonuniform, but a rule was 
followed that no interval be more than twice as 
wide as its neighbor. 

The source was represented by a top bound- 
ary source in the first six downward angles. The 
source varied with energy group according to the 
measured source spectrum and with radial inter- 
val according to the measured beam profile. 

In the DUCT calculations, the duct was re- 
presented by a cylinder 2. 54 cm in radius com- 
pletly penetrating the shield with its axis coincid- 
ing with the axis -of- symmetry of the shield. The 
rest of the geometry of the shield was input to the 
GEOM package, with the requirement that the duct 
be re-input as the first quadratic surface for 
GEOM. 

DUCT calculations were run for three de- 
tector positions with shield, sodium, or void in 
the duct. In each case the objective is to calculate 
the flux crossing the duct wall or colliding within 
the duct and reaching the detector. 

To "know" the unperturbed flux at the de- 
tector location requires that it be available from 
the DOT calculation. For detectors at or in the 
shield, simple interpolation of the DOT fluxes will 
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suffice. For detectors near the shield it is there- 
fore recommended that the DOT calculation con- 
tain void intervals which themselves contain the 
proposed detector sites. For detectors far from 
the shield, surface integration of the DOT fluxes 
and their projection to the detector, such as is 
done in DASH is required to determine the un- 
perturbed flux. 

3. 1 COMPARISON WITH MONTE CARLO 

CALCULATION 

A MORSE (ref. 4) Monte Carlo calculation 
was performed for the experimental configuration 
with sodium in the duct. Detector No. 1 was re- 
presented by a point one inch from the shield on 
the duct axis. The scattered flux was calculated 
by last-flight estimation from each collision site. 
The source was represented by a circular beam 
six inches in radius with its center on the duct 
axis. The general 05R geometry package (GEO M) 
was used. 

The final calculation involved 84 batches 
of 100 neutrons each. Each neutron started from 
some point in the 6-in. radius circle on the top 
surface of the shield. Its distance from the cen- 
ter of the circle was selected uniformly from 
0 in. to 6 in. and its weight adjusted by 

WATE = 2*WATE*RADIUS/BEAM, 

where RADIUS is its distance and BEAM = 6 in. , 
to account for the fact that the circular area was 
not sampled uniformly. Its direction from the 
center of the circle was selected randomly. Its 
starting direction cosines were (0, 0, 1), that is, 
normal incidence. 

The initial energy group of the neutron 
was selected from the 15 groups not according to 
the natural source energy distribution but from a 
biased distribution which favored the higher en- 
ergy groups; MORSE adjusted its weight to com- 
pensate for the source biasing. 

The path length of each flight was selected 
not from the natural exponential distribution but 
from a transformed exponential distribution which 
encouraged longer paths near the the favored 
downward direction and shorter paths near the 
upward direction. The parameter PATH = 0. 5 
was used; MORSE adjusts the weight to account 
for the path length biasing. 

At the end of each flight path is a collision 
site; since absorption is not allowed, the neutron's 
weight is suitably multiplied by a nonabsorption 
probability. A last-flight estimate is then made, 
contributing 

P(0)*WATE*EXP(-ARG)/R 2 

to the scalar scattered flux at the detector for 
every energy group into which the neutron could 
possibly scatter from its present group. Then 
the neutron is allowed to launch its next flight in 
a newly determined direction with path length 
selected from the transformed exponential distri- 
bution. 

At the end of the 84 batches, the scattered 


flux and its fractional standard deviation were 
calculated for each energy group, the scattered 
flux was added to the direct flux, and the results 
divided by the group width in MeV. The final units 
of the flux were (neutrons /cm 2 -MeV) /(Incident 
neutron-cm' 2 ). Eighty thousand, four-hundred 
eighty- two scattering events contributed to the flux 
estimates. The fractional standard deviation of the 
scattered flux ranged from 5 to 20 percent. 

A comparison with the DUCT calculation is 
shown in fig. 1. In general, the agreement varies 
between 10 to 20 percent, with the DUCT calcula- 
tions consistently higher, except for the lowest 
two energy groups. 



Fig. I. Comparison between DUCT and MORSE calculated 
axial responses 


4. EXPERIMENTAL PROGRAM 

The experiments were designed to be as 
sensitive as possible to streaming but "clean" 
enough to provide a good check on the analysis. 

4. 1 NEUTRON SOURCE 

The experimental layout for the experi- 
ments is shown in fig. 2. Photo-neutrons were 
produced in the tantalum-beryllium converter bom- 
barded by 50- MeV electrons from the Gulf Radia- 
tion Technology linear accelerator. This conver- 
ter has been used in several previous programs 
(ref. 5) at Rad Tech and provides the source spec- 
trum shown in fig. 3. This is the spectrum mea- 
sured by time of flight over 50 meters; 2. 54 cm of 
lead and 3.144 cm of uranium were used as filters 
at the 16-meter station to reduce the gamma flash. 
This setup is for the lithium -hydride- sodium duct 
experiments. A similar setup was used for the 
BATH-NERVA shield streaming experiments with 
the exception that the measurements were made at 
16. 7 meters. Neutrons produced in the tantalum- 
beryllium target were then collimated to provide a 
circular plane beam for normal incidence on the 
shield assembly. 
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Fig. 2. Experimental layout for the lithium hydride duct 
shield experiment 


Fig. 4. Lithium Hydride Duct Shield Assembly, 


Fig. 3. Measured source spectrum for the lithium hydride 
duct shield experiment 


4. 2 SHIELD ASSEMBLIES 


Two shield assemblies were studied. The 
first assembly was cylindrical lithium hydride 
assembly 50. 8 cm in diameter by 25. 4 cm thick 
with a 5. 08 cm diameter axial duct. The assembly 
in situ at the end of the 50 -meter flight path is 
shown in fig. 4. The duct was filled with either 
sodium or lithium hydride. 

The BATH assembly consists of nineteen 
two-in. hexagonal by 11-in. long pieces, placed in 
an hexagonal matrix approximately 10 inches dia- 
meter by 11.1 inches long. To minimize the effect 
of room-return neutrons, the matrix was enclosed 
in a 10-in. diameter cylindrical container enclosed 
by a 5-in. thick shield. The shield consisted of a 
mixture of 50% by weight of lithium carbonate 
(Li 2 C 03 ) and polyethylene (CH 2 ). The inner BATH 
assembly is shown in fig. 5; the assembled BATH 
and side shield in situ at the end of the 16-meter, 
flight path is shown in fig. 6. 


Figure 6. 


Inner- BATH Shield Assembly. 

',r- 


BATH Shield Assembly Enclosed 
in Outer Shield 




4.3 


DATA ACQUISITION 


During a previous Rad Tech program (ref. 
5), a new experimental method was developed 
for studying simultaneously fast-neutron trans- 
port and secondary gamma- ray production and 
transport in shields. With this method, an elec- 
tron accelerator is used to produce intense 50- 
nsec pulses of photoneutrons. These neutrons 
stream down an evacuated 52-meter long drift 
tube to strike a shield assembly. An NE-213 de- 
tector is positioned immediately behind the shield 
assembly to detect both fast neutrons and second- 
ary gamma rays produced in a shield assemblyby 
(n, xy) and (n, xy) reactions. 

A simplified block diagram of the three- 
parameter data acquisition system used with the 
NE-213 detector is shown in fig. 7. Three para- 
meters, t, h, and x, are recorded for each detec- 
tion event. The time parameter t is the time in- 
terval between the Linac pulse and the detection 
event; the pulse-height parameter h is proportion- 
al to the light produced in the NE-213 scintillator 
by neutron and gamma- ray interactions; and the 
pulse- shape discrimination parameter x identifies 
detection events as due to neutrons or gamma 
rays. The key features of each parameter are 
discussed next. 



Fig. 7. Simplified diagram of the three -parameter data 
acquisition system 


4. 3. 1 Time Parameter 

The time parameter is related to the 
energies of neutrons incident on the shield assem- 
bly rather than the energies of neutrons leaking 
from the shield. This is true because the shield 
assembly is located at the end of the flight path 
near the detector rather than at the beginning 
near the pulsed neutron source as it would be for 
a conventional time -of -flight shielding experi- 
ment. The present geometry arrangement is also 
used for (n, xy) and (n, y) reaction studies (ref. 6). 

The time scale may be naturally divided 
into two regions. For short times (t < 4 fj, sec), 
detection events due to: (1) neutrons with energies 
above the detector threshold at 0. 9 MeV (t = 4. 0 
/jsec), and (2) gamma rays produced by (n, xy) 
reactions with threshold at 1. 3 MeV (t = 3. 3 /nsec) 
are recorded. These events may be analyzed and 
interpreted in either the time domain or incident 
neutron energy domain. The results are com- 
pletely equivalent within the abovementioned 
appr oximation. 

However, for long times (t > 4 yusec), neu- 
trons are not detected and only gamma rays pro- 
duced by (n, y) reactions are detected. Because 
these reactions may be due to uncollided neutrons, 
slowing-down neutrons, and/or thermalized neu- 
trons, time and incident neutron energy cannot be 
uniquely related. Thus, these data are analyzed 
only in the long time domain. In summary, the 
data which are obtainable are categorized as due 
to: 

1. Short time or fast neutrons 

(E > 0. 8 MeV, t < 4. 0 /usee) 

2. Short time (n, xy) or inelastic gamma 

rays (t < 3. 3 jusec) 

3. Long time (n, y) or capture gamma rays 

(t > 3. 3 jusec). 

4. 4 UNFOLDING 

The mathematical procedure for unfolding 
i6 carried out in the FERDOR code. The method 
is discussed in refs. 7 and 8. The results pre- 
sented here are only the output from the FERDOR 
code. Each unfolded spectrum is not uniquely 
defined but instead gives a band of values within 
which the actual spectrum should lie. 

The calculation of this band by FERDOR 
includes two contributions, one due to the statisti- 
cal errors in the input pulse-height data, and the 
other to the suboptimal calculation of the "inverse" 
of the response matrix, because, in general, such 
errors are: (1) difficult to assess due to the indi- 

rect manner in which response matrices are pro- 
duced by interpolation from calculation, and a 
measurement for example, and (2) difficult to pro- 
pagate through the unfolding procedure. In our 
case, the response matrix was obtained entirely 
from a Monte Carlo 06R calculation performed at 
the Oak Ridge National Laboratory. 
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An excellent check on the methods des- 
cribed in the previous section is provided by the 
time -of -flight spectrum. Note that this spectrum 
is always measured simutlaneously by the three- 
parameter data acquisition system. The mea- 
sured source spectrum (no shield assembly) ob- 
tained by time of flight should agree with that 
obtained by unfolding the proton pulse height spec- 
tra. A comparison is shown in fig. 8; the agree- 
ment is very good, which lends confidence to the 
unfolding procedure used here. 



Fig. 8. Comparison between the FERDOR unfolded and 
time-of-flight measured source spectra 

5. COMPARISON OF DUCT CALCULATION 
AND EXPERIMENT 

Comparisons between the experimental data 
obtained with the lithium hydride assembly and 
calculations performed with the DUCT code are 
shown in figs. 9 . 10, and 11. The calculations 
yielded neutron flux per MeV per source neutron, 
the measurements yielded neutron flux per MeV 
per source monitor. To compare the two on a 
source neutron basis, the measured results were 
multiplied by a source normalization factor. This 
factor was obtained by integrating the calculational 
source spectrum and the measured source spec- 
trum over the same energy limits; namely, 2. 0 
MeV to 5. 5 MeV, and then dividing the calculation- 
al source integral by the experimental source 
integral, that is 


IV* 5 1 

- / <p(E)dE 

Jz.o 

DOT 

V 5 ' 5 


/ tp(E)dE 


.2. 0 

EXP 



Fig. 9. Comparison between measured and DUCT calculated axial 
detector responses for a detector located at the duct exit 



Fig. 10. Comparison between measured and DUCT calculated axial 
detector responses for a detector located at the duct exit 



Fig. II. Comparison between measured and DUCT calculated axial 

detector responses for a detector located 15. 24 cm from the duct exit 
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This was done for the source spectrum only so ' 
that no further normalization is applied when 
comparing measurement with theory. 

The comparison between the measured and 
axial detector responses is shown in fig. 9 for 
the sodium and lithium hydride filled ducts. 
Figure 10 compares the measured and calculated 
fluxes for void and lithium hydride filled ducts. 
Figure 11 compares the calculated and measured 
axial detector respones for a collimated detec- 
tor located 12. 7 cm from the duct exit. The 
final results for the BATH experiment-calcula- 
tion comparisons are not complete at this 
writing. 

6 . CONCLUSIONS 

In view of the overall agreement between 
measurement, DUCT and MORSE, we feel the 
DUCT code provides a good engineering code to 
claculate the effect of coolant type holes in 
shields. We have tested the method to its ex- 
treme limit for streaming. 

The experiments were limited to neutron 
energies above 1 MeV; however, as noted earlier, 
the agreement with MORSE was very good at the 
low energies (E <10 keV) where the fluxes tend to 


be isotropic lending confidence to the DUCT cal- 
culations at these energies. It is important to 
note that the MORSE calculation required approx- 
imately 30 to 45 minutes to achieve a 10 to 20% 
variance on the scattered flux estimates for a 
single detector position. The DUCT calculations 
require approximately four minutes per detector 
with an additional 40 to 60 minutes for the DOT 
calculation. However, if one were making an 
engineering type analysis of a shield where many 
duct-detector configurations are to be studied, 
the DOT-DUCT running time would be consider- 
ably shorter than the MORSE running time, as 
seen in the following table: 

MORSE VS DOT /DUCT RUNNING 
TIMES PER DETECTOR (min) 

Number of Ducts 12 4 6 

MORSE 40 80 160 240 

DOT/DUCT 64 68 72 80 

We recommend that a void region be in- 
cluded in the DOT calculation for detectors near 
the shield/void region to minimize errors in ex- 
trapolating the DOT fluxes. For detectors quite 
distant from the shield, we recommend the use of 
DASH, which we plan to incorporate in the pro- 
duction version of DUCT. 
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